/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package provider.impl.genomemap.lc;

import commons.core.LikelihoodComputer;
import genomemap.data.TAB08Data;
import genomemap.data.TAB08Data.Genotype;
import org.apache.commons.lang.ArrayUtils;

import java.math.BigDecimal;

/**
 * @author Susanta Tewari
 * @since Aug 21, 2011
 */
class TAB08LCSFA implements LikelihoodComputer<TAB08Data> {

    /**
     * Method description
     *
     *
     * @param s
     *
     * @return
     *
     * @throws InterruptedException
     */
    @Override
    public BigDecimal computeLogLHood(final TAB08Data s) throws InterruptedException {

        final byte[][] f = new byte[s.genotypeCount()][];
        final int[] n    = new int[s.genotypeCount()];
        int counter      = 0;

        for (final Genotype genotype : s) {

            f[counter] = genotype.getScore();
            n[counter] = genotype.getFreq();

            counter++;
        }

        final double[] c =
            ArrayUtils.toPrimitive(s.getModel().getExchangeProbs().toArray(new Double[0]));
        final double[] p = implementHMM(f, c);
        double logLHood  = 0;

        for (int i = 0; i < n.length; i++) {
            logLHood += n[i] * Math.log(p[i]);
        }

        return new BigDecimal(logLHood);
    }

    double[] implementHMM(final byte[][] f, final double[] c) throws InterruptedException {

        final double[] lhood = new double[f.length];

        /* x is crossover as defined in model description */
        final ExchangeConfig[] x = new ExchangeConfig[c.length];

        for (int i = 0; i < x.length; i++) {
            x[i] = new ExchangeConfig();
        }

        do {

            for (int i = 0; i < f.length; i++) {
                lhood[i] += getConditionalProb(f[i], x) * getMarginalProb(c, x);
            }

            if (Thread.currentThread().isInterrupted()) {
                throw new InterruptedException();
            }

        } while (getCrossoverPattern(x));

        return lhood;
    }

    double getConditionalProb(final byte[] f, final ExchangeConfig[] x) {

        byte[] prev_strand        = { 1, 1, 0, 0 };
        final byte[] strand       = prev_strand.clone();
        final byte[] strand_match = new byte[4];

        if (f[0] == 1) {
            strand_match[0] = strand_match[1] = 1;
        } else if (f[0] == 0) {
            strand_match[2] = strand_match[3] = 1;
        }

        for (int i = 0; i < x.length; i++) {

            modifyStrand(strand, x[i].getI());
            modifyStrand(strand, x[i].getJ());

            for (int j = 0; j < 4; j++) {

                if ((prev_strand[j] != f[i]) || (strand[j] != f[i + 1])) {
                    strand_match[j] = 0;
                }
            }

            if (sum(strand_match) == 0) {
                return 0;
            }

            prev_strand = strand.clone();
        }

        return (double) sum(strand_match) / 4.0;
    }

    /**
     * @param strand genotype strand
     * @param i a single exchange.i or j in ExchangeConfig for any interval
     */
    protected void modifyStrand(final byte[] strand, final byte i) {

        switch (i) {

            case 1 :
                swap(strand, 0, 2);

                break;

            case 2 :
                swap(strand, 1, 2);

                break;

            case 3 :
                swap(strand, 1, 3);

                break;

            case 4 :
                swap(strand, 0, 3);
        }
    }

    double getMarginalProb(final double[] c, final ExchangeConfig[] x) {

        double prod = 1.0;

        for (int i = 0; i < x.length; i++) {

            if (x[i].getI() == 0) {
                prod *= 1 - c[i];
            } else {
                prod *= c[i] / 4;
            }

            if (x[i].getJ() == 0) {
                prod *= 1 - c[i];
            } else {
                prod *= c[i] / 4;
            }
        }

        return prod;
    }

    boolean getCrossoverPattern(final ExchangeConfig... x) {

        for (int i = x.length; --i >= 0; ) {

            if (x[i].nextPattern()) {
                return true;
            } else {
                x[i] = new ExchangeConfig();
            }
        }

        return false;
    }

    /**
     * Method description
     *
     *
     * @param array
     * @param i
     * @param j
     */
    private void swap(final byte[] array, final int i, final int j) {

        final byte temp = array[i];

        array[i] = array[j];
        array[j] = temp;
    }

    static byte sum(final byte[] array) {

        byte result = 0;

        for (final byte val : array) {
            result += val;
        }

        return result;
    }
}
